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Abstract 

Large-Eddy  Simulation  of  transcritical  turbulent  mixing  is  investigated  in  this 
report.  This  analysis  is  motivated  by  rocket  propulsion  applications  where  cryo¬ 
genic  oxygen  is  injected  at  low  temperature,  below  the  critical  value,  in  a  high 
pressure  environment  where  the  pressure  is  above  critical.  This  situation  is  here 
envisaged  by  analyzing  the  behavior  of  a  nitrogen  jet  delivered  by  a  shear  coaxial 
injector  comprising  an  inner  stream  injected  at  temperatures  close  to  the  critical 
temperature  and  a  coaxial  flow  at  a  higher  temperature.  Real  gas  Large  Eddy  Sim¬ 
ulations  are  carried  out  for  selected  experiments  previously  performed  at  AFRL, 
Edwards  (USA).  Two  geometrical  configurations  associated  with  two  sets  of  mass 
flow  rates  have  been  studied  at  AFRL  (LarThin  and  SarThick).  They  lead  to  dif¬ 
ferent  jet  dynamics  and  subsequently  different  responses  to  acoustic  modulation. 
Calculations  retrieve  what  is  observed  experimentally  both  in  the  absence  of  an 
external  modulation  and  under  transverse  acoustic  field.  It  is  shown  in  particular 
that  in  the  LarThin  case  the  acoustic  forcing  produces  large  scale  transverse  oscil¬ 
lations  of  the  inner  stream  and  promotes  the  formation  of  counter  rotating  vortices 
in  the  wake  of  the  high  density  jet.  These  coherent  structures  stretch  the  inner 
jet,  increasing  the  mass  exchange  surface  between  the  two  streams  and  the  mixing 
efficiency.  In  the  SarThick  situation,  a  large  scale  back  flow  is  established  around 
the  inner  jet,  which  limit  its  downstream  penetration  promoting  recirculation  of 
cold  fluid  behind  the  thick  lip.  The  acoustic  modulation  has  a  limited  effect  on  the 
jet  geometry  in  this  case,  probably  because  of  the  presence  of  this  dynamic  back 
flow  region.  The  inner  jet  length  is  nevertheless  reduced  in  a  significant  way. 

Keywords  ;  Supercritical  injection ,  shear  coaxial  injectors,  transverse  acoustic  modula¬ 
tion. 


1  Introduction 

Modern  liquid  rocket  engines,  which  use  cryogenic  propellants,  commonly  operate  under 
transcritical  conditions,  a  situation  where  the  chamber  pressure  is  above  the  critical  point 
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and  the  injection  temperature  is  below  the  critical  value  for  at  least  one  of  the  propellants. 
During  the  design  phase  of  these  engines,  combustion  instabilities  (Cl)  may  occur,  which 
can  stall  or  severely  delay  their  development.  CIs  result  from  a  complex  coupling  between 
acoustic  modes  of  the  chamber  and  the  flow,  mixing  and  heat  release  processes.  The 
latter  is  usually  structured  by  the  inner  and  outer  streams  originating  from  a  shear 
coaxial  element  (most  common  injector  design  for  cryogenic  engines),  interacting  with 
each  other  and  with  their  surroundings.  The  dynamics  and  mixing  of  coaxial  transcritical 
jets  submitted  to  acoustic  modulations  is  then  of  interest  for  the  study  of  combustion 
instabilities.  Moreover,  studying  the  transient  and  unsteady  phenomena  occurring  in 
such  configurations  is  of  fundamental  importance  for  the  understanding  and  prediction 
of  cryogenic  flow  dynamics.  This  is  here  investigated  by  analyzing  the  behavior  of  a 
nitrogen  shear  coaxial  jet  comprising  an  inner  stream  injected  at  temperatures  close 
to  the  critical  temperature  and  a  coaxial  flow  at  a  higher  temperature.  Experiments  are 
carried  out  both  in  the  absence  of  external  modulation  and  by  imposing  a  large  amplitude 
transverse  acoustic  field.  Real  gas  Large  Eddy  Simulations  (LES)  are  performed  for 
selected  experiments.  The  combination  of  experiments  and  calculations  is  used  to  evaluate 
effects  of  the  operating  parameters  on  the  coaxial  jet  response. 

The  interaction  of  transverse  acoustic  waves  with  transcritical  flows  and  its  influence  on 
mixing  efficiency  has  been  considered  in  previous  studies.  The  impact  of  acoustics  on  a 
non-reacting  shear  coaxial  jet  has  been  investigated  experimentally  by  Leyva,  Rodriguez 
et  al.  [15,  14,  4],  They  found  that  the  injector  geometry  has  an  effect  on  the  sensitivity  of 
the  coaxial  flows  to  acoustic  modulation.  This  study  has  been  continued  recently,  through 
a  grant  from  EOARD,  by  considering  the  combination  of  experiments  and  simulations,  in 
a  collaboration  between  AFRL  and  EM2C  [16].  This  collaboration  was  quite  productive 
in  terms  of  shedding  light  on  key  experimental  observations  made  at  AFLR  by  using 
state  of  the  art  LES  tools  available  at  EM2C  (see  [17,  18]  for  an  example  of  applications). 
Simulations  were  carried  out  with  the  AVBP-RG  software  jointly  developed  by  Cerfacs 
and  EM2C.  By  analyzing  the  spectral  content  of  the  LES  flow  field  results,  it  was  found 
that  in  this  particular  configuration  the  transcritical  inner  jet  was  relatively  insensitive 
to  the  acoustic  modulation,  while  the  gaseous  annular  jet  was  strongly  perturbed  and 
generated  large  scale  vortices.  It  was  suggested  that  the  relative  insensitivity  of  the  inner 
flow  is  linked  to  the  fact  that  the  frequency  was  above  the  range  of  preferred  frequencies 
of  the  central  jet.  Key  natural  hydrodynamic  instabilities  were  also  determined  with 
LES,  which  is  of  paramount  importance  to  calculate  which  frequencies  can  interact  with 
the  combustion  process. 

The  general  goal  is  to  continue  these  investigations  on  the  mixing  of  coaxial  transcriti¬ 
cal  jets  submitted  to  acoustic  modulations.  In  the  previous  collaborative  study  between 
AFRL  and  EM2C,  a  single  injector  geometry,  with  a  thick  inner  jet  lip  separating  the 
2  streams,  was  considered.  One  particular  objective  of  the  present  research  is  to  assess 
the  influence  of  geometrical  parameters  on  the  jet  dynamics  and  sensitivity  to  acoustic 
modulation  by  using  experimental  data  from  AFRL  on  an  alternate  geometry,  a  thin  lip 
injector  [15].  Key  experimental  cases  will  be  simulated  using  LES.  It  has  been  experi¬ 
mentally  observed  at  AFRL  that  the  flow  in  this  alternative  geometry,  when  exposed  to 
acoustics,  shows  more  visible  effects  such  as  a  large  shortening  of  the  inner  jet  “intact” 
length.  LES  can  help  identify  the  mechanisms  which  define  this  dynamical  behavior. 
Both  experiments  and  numerical  simulations  are  used  in  the  present  study.  Experiments 
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are  carried  out  at  AFRL  while  numerical  simulations  are  performed  at  the  EM2C  labora¬ 
tory.  Experiments  provide  high  speed  shadowgraphs  of  the  flow  and  average  temperature 
profiles.  From  these  data,  simulations  can  be  validated  and  results  of  calculations  can 
then  be  used  to  augment  the  experimental  measurements  and  visualizations  in  order  to 
study  the  flow  dynamics  and  identify  the  subtending  processes. 


2  Balance  equations  and  models 

This  section  provides  a  synthesis  of  the  main  equations  and  submodels  used  in  the  present 
simulations. 


2.1  Governing  equations 

A  single  species  flow  configuration  is  investigated  in  the  present  investigation.  In  this 
situation,  the  vector  of  conservative  variables  for  a  compressible  flow  is  w  =  (p,  pu,  pE)T , 
where  p  is  the  density,  u  the  velocity  vector  and  E  the  total  energy  (E  is  the  sum  of 
internal  energy  es  and  kinetic  energy  ep:  =  1/2^  ) .  The  large-eddy  simulations  are 

carried  out  by  integrating  the  mass  weighted  spatially  filtered  Navier-Stokes  equations. 
In  the  following  balance  equations  (for  a  given  variable  0,  (p  designates  the  spatially 
filtered  variable,  while  <j>  is  the  mass-weighted  (Favre)  spatially  filtered  variable)  : 


dpE 

~df 


dp 

dt 

d~pui 

dt 

dpEUj 

dxj 


dpUj 
+  dxj 

=  0 

(1) 

dpujUi 

dp 

df-j 

dr- 

ij 

(2) 

dxj 

dxi 

dxj 

dxj 

dpUj 

dq) 

dUiT~j 

(3) 

+  9xj 

dxj 

dXj 

dxj 

the  spatial  coordinates  vector  and  t  is  the  time. 

The 

In  Eqs.  1-3,  p  is  the  pressure, 
laminar  viscous  stress  tensor  r  and  heat  flux  vector  q  are  expressed  as  linear  functions 
of  the  strain  rates  and  temperature  gradient  respectively.  Laminar  viscosity  and  heat 
conductivity  coefficients  are  determined  with  the  Chung  et  al.  method  [2] .  For  the  single 
species  case  investigated  in  the  present  article,  no  additional  diffusion  terms  are  needed 
since  the  Soret  and  Dufour  effects  naturally  vanish.  The  sub  grid  scale  tensor  rk  rely  on 
the  WALE  model  [7]  and  the  subgrid  scale  heat  flux  is  modeled  using  a  standard  gradient 
transport  assumption  : 


**  =  <4) 

The  subgrid  scale  conductivity  is  cast  in  the  form  At  =  cppt/Pvt  where  pt  =  put  and  Prt 
respectively  designate  the  subgrid  dynamic  viscosity  and  turbulent  Prandtl  number.  As 
in  low  pressure  simulations  a  constant  value  Prt  =  0.7  is  used  in  the  present  calculations. 
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2.2  Real  gas  thermodynamics  and  equation  of  state 

Real  gas  thermodynamics  for  high  pressure  flow  calculations  is  envisaged  in  many  recent 
studies  where  accurate  descriptions  of  the  fluid  state  are  derived  from  various  types  of 
equations  [1,  21].  The  present  simulations  rely  on  the  Peng- Robinson  equation  of  state 
(EOS)  [10]  (Eq.  5)  which  is  used  in  the  analysis  to  model  the  departure  from  an  ideal-gas 
behavior.  The  fluid  state  is  expressed  in  the  form  : 


prT  p2a(T) 

1  —  pb  1  +  2  pb  —  p2b2 


(5) 


where  T  is  the  temperature,  r  =  R/W  with  R  the  perfect-gas  constant  and  W  the  molar 
mass.  The  Peng- Robinson  coefficients  a(T )  and  b  for  a  single-species  fluid  are  [12]: 


a 

b 


0.457236^^-  (l  +  c(l 

Pc  V 
vT 

0.077796— 

Pc 


(6) 

(7) 


with 


c  =  0.37464  +  1.54226a;  —  0.26992a;2  (8) 

where  Tc  =  126.2  K  and  pc  =  3.396  MPa  are  the  critical  temperature  and  pressure  of 
nitrogen,  u  =  0.0372  represents  its  acentric  factor  and  Tr  =  T/Tc  is  the  reduced  temper¬ 
ature.  This  equation  offers  a  good  trade-off  between  computational  cost  and  precision. 
The  pressure  dependence  of  thermodynamic  functions  and  coefficients  is  subsequently 
derived  from  the  EOS  [5,  8,  21]. 


3  Computational  setup 

3.1  Numerical  aspects 

The  mass  weighted  spatially  filtered  compressible  Navier-Stokes  equations  in  combina¬ 
tion  with  the  subgrid  scale  models  described  previously  are  integrated  in  the  AVBP  flow 
solver  [19,  6].  The  present  integration  method  relies  on  a  Taylor-Galerkin  weighted  resid¬ 
ual  central  distribution  scheme,  called  TTGC.  This  scheme  is  third-order  in  space  and 
third-order  in  time  [3]. 

The  real  gas  equation  of  state,  thermodynamics  and  transport  coefficients  have  been 
implemented  in  AVBP.  The  convective  fluxes  Jacobian  matrices,  used  by  the  scheme, 
are  expressed  in  terms  of  real  gas  thermodynamics  to  preserve  the  overall  consistency 
of  the  code.  A  fully  consistent  treatment  of  associated  boundary  conditions  is  based 
on  the  characteristic  wave  decomposition  method  NSCBC  [11,  6]  and  includes  specific 
expressions  provided  for  real  gas  thermodynamics  in  [9] .  The  treatment  of  steep  gradients 
relies  on  artificial  viscosity  and  is  activated  when  unresolved  gradients  are  detected  by 
a  density-based  sensor.  More  details  about  the  sensor  and  the  operators  used  to  trigger 
artificial  viscosity  addition  can  be  found  in  [16]. 
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3.2  Computational  domain 

The  mixing  chamber  used  in  the  simulation  is  shown  in  Fig.  la.  Its  dimensions  corre¬ 
sponds  to  the  experimental  ones,  with  a  chamber  length  of  59.4  mm,  and  a  rectangular 
section  of  12.7  x  76  mm.  In  the  experiment,  in  order  to  maintain  the  pressure  in  the 
chamber,  the  chamber  is  placed  inside  a  larger  pressurization  system.  The  two  chambers 
are  connected  through  an  upstream  perforated  plate  and  through  the  outlet  pipe.  The 
velocity  through  the  perforated  plate  is  thus  not  imposed,  but  results  from  entrainment 
by  the  coaxial  jet.  This  geometry  can  be  qualitatively  accounted  for  in  the  simulation  by 
placing  an  upstream  large  volume  as  shown  in  Fig.  lb. 

In  order  to  assess  the  influence  of  the  perforated  plate  holes,  preliminary  simulations  were 
carried  out  with  and  without  the  upstream  reservoir.  When  the  the  upstream  reservoir 
is  not  used,  an  inlet  velocity  is  imposed  at  the  perforated  walls.  Two  values  for  the  inlet 
velocity  are  compared:  1  m  s-1  and  0.1  m  s-1  to  test  the  sensitivity  to  this  parameter. 
Results  are  plotted  in  Fig.  2  for  a  set  of  radial  profiles.  This  shows  that  the  profiles 
are  close  to  each  other.  Thus  an  injection  velocity  can  be  imposed  at  the  hole  exhaust 
section,  which  simplifies  the  computational  domain  and  lowered  the  mesh  requirement. 
A  value  of  1  m  s-1  is  chosen  in  the  calculations. 


Figure  1:  Computational  domain  (a)  Mixing  chamber  (b)  Domain  with  the  upstream 
reservoir  connected  to  the  perforated  plate.  Preliminary  simulations  show  that  the  up¬ 
stream  reservoir  can  be  neglected  for  the  simulations,  and  a  low  velocity  injection  can  be 
imposed  at  the  perforated  plane  inlets  instead. 
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Figure  2:  Influence  of  perforated  walls  boundary  conditions.  Imposed  inlet  velocity:  o 
0.1  m/s  —  1  m/s  ;  with  an  upstream  reservoir  . 


3.3  Injector  geometries  and  mesh  definitions 

Two  injector  geometries  are  considered  in  this  study.  They  differ  by  the  thickness  of  the 
inner  lip  and  the  surface  area  ratio  between  the  outer  and  inner  injection  plane.  They  are 
referred  to  as  LarThin  (Large  aspect  ratio,  Thin  lip)  and  SarThick  (Small  aspect  ratio, 
Thick  lip).  Their  characteristics  are  given  in  Tab.  1. 


Case 

Si  [mm] 

Se/Si 

[mm] 

he  [mm] 

LarThin 

SarThick 

0.095 

1.13 

10.6 

2.9 

0.7 

1.47 

0.775 

0.37 

Table  1:  Injector  geometrical  parameters.  8i  is  the  inner  lip  thickness,  Se  the  outer 
injector  section  surface,  S,  the  inner  post  section  surface,  d,  the  inner  injector  diameter 
and  he  the  outer  injector  thickness. 


Figure  3:  Longitudinal  slice  of  the  meshes  (near  injector  region)  (a)  LarThin  (b)  SarThick. 
SeeTab.  2  for  details 


Longitudinal  slices  of  the  meshes  for  cases  LarThin  and  SarThick  are  shown  on  Fig.  3 
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and  4.  They  are  both  fine  near  the  injector  with  a  smooth  coarsening  toward  exit.  For 
SarThick  simulations,  the  mesh  is  refined  in  the  annular  injector  in  order  to  ensure  a 
proper  turbulence  resolution.  The  grid  characteristics  are  detailed  in  Tab.  2.  Both 
meshes  comprise  approximately  20  million  cells,  with  30  cells  along  the  inner  and  outer 
injector  characteristic  length  for  case  LarThin,  and  20  and  15  cells  in  the  inner  and  outer 
injectors  for  case  SarThick.  A  grid  sensitivity  study  is  proposed  in  Sec.  5. 


Figure  4:  Longitudinal  slice  of  the  mesh  (a)  LarThin  (b)  SarThick  (Tab.  2) 


Table  2:  Characteristics  of  the  mesh  used  for  the  simulations.  At  is  the  time  step, 
rc  =  25 dinner / Uinner  the  convective  time. 

3.4  Injection  characteristics  and  boundary  conditions 

These  simulations  consider  pure  nitrogen  injected  in  a  chamber  at  a  pressure  of  3.56  MPa 
for  both  cases.  The  chamber  temperature  is  223  K  for  case  LarThin  and  205  K  for  case 
SarThick.  Injection  conditions  for  the  simulated  cases  are  given  in  Tab.  3  and  4.  The 
two  cases  have  comparable  density  ratio  between  the  outer  and  inner  injection  channels 
(7.7  for  LarThin  and  9.0  for  SarThick).  There  is  a  slight  departure  between  the  two  cases 
considering  the  velocity  ratio  (6.5  for  LarThin  and  9.0  for  SarThick),  which  leads  to  a 
higher  momentum  flux  ratio  for  case  SarThick  than  for  LarThin  (9.2  vs  5.0). 

All  the  solid  boundaries  in  the  simulation  are  treated  as  no  slip  adiabatic  walls.  Pressure 
is  imposed  at  the  outlet  using  a  partially  non-reflecting  characteristic  boundary  condition. 
Flat  profiles  are  injected  at  the  perforated  plate.  Fully  turbulent  profiles  are  assumed 
at  injector  inlets  with  anisotropic  turbulent  fluctuations  superimposed  to  the  mean  flow 
profile  at  the  injector  inlet  (Fig.  5)  [20]. 
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Case 

T ?  [K] 

u* 

m  s  1 

pi  [kg  m  3 

Te  [K] 

ue 

ms1 

Pe  [kg  m  3 

LarThin 

SarThick 

122 

120 

2.5 

2.1 

531 

556 

203 

203 

16.2 

18.9 

63 

63 

Table  3:  Injection  characteristics  of  the  simulated  cases,  pi,  pe,  Ui,  ue,  Tj  and  Te  are  the 
inner  and  outer  jet  densities,  velocities  and  temperature. 


Case 

J  =  (peu2e)/{piu2i) 

M  =  (peue)  /  (piUi) 

S  pi/  pe 

U  =  ue/ui 

LarThin 

5.0 

0.77 

7.7 

6.5 

SarThick 

9.2 

0.72 

9.0 

9.0 

Table  4:  Injection  characteristics  of  the  simulated  cases.  J  is  the  momentum  flux  ratio 
between  the  outer  and  the  inner  jets,  M  the  mass  flux  ratio,  S  the  density  ratio  and  U 
is  the  velocity  ratio,  pi,  pe,  and  ue  are  the  inner  and  outer  jet  densities  and  velocities. 


Figure  5:  Radial  velocity  profiles  at  injection  for  both  inlet,  (a)  Axial  mean  velocity  u, 
normalized  by  its  centerline  value  uc,  (b)  axial  ,  azimuthal  and  radial  rms  velocities  (  u'x, 
u'i  and  u'r,  respectively). 
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4  Grid  sensitivity  study 

To  estimate  the  grid  sensitivity  calculations  are  carried  out  for  the  two  geometries  on  three 
different  meshes  referred  as  “coarse” ,  “intermediate”  and  “fine” .  Their  characteristics  are 
detailed  in  Tab.  5.  They  all  share  the  same  cell  distribution,  shown  in  Fig.  4  and  only 
differ  by  the  smallest  characteristic  cell  size  Axmin  present  in  the  grid.  Comparisons 
between  results  obtained  with  the  three  meshes  indicates  that  the  intermediate  mesh 
is  adequate.  This  mesh  is  used  in  the  final  simulations  in  this  document.  This  study 
requires  large  CPU  resources  and  it  is  important  to  sure  make  sure  that  the  calculations 
are  accurate  but  with  a  minimal  expenditure  of  CPU  time. 


Case 

Nb  of  cells 

At  [s 

CPU  for  1  rc  [h] 

Axmin  [mm] 

LarThin  coarse 
LarThin  intermediate 
LarThin  fine 

6  237  282 
20  944  365 
77  301  292 

2.9  10"8 

1.9  10“8 
1.05  10“8 

2  000 

10  200 

68  000 

0.0375 

0.025 

0.017 

SarThick  coarse 
SarThick  intermediate 
SarThick  fine 

8  180  162 
23  664  745 
64  763  744 

2.8  10“8 
1.9  10"8 
1.4  10"8 

11  600 

38  750 

151  000 

0.035 

0.0245 

0.0175 

Table  5:  Grid  characteristics  of  the  simulated  cases.  At  is  the  time  step,  rc  =  25d, /ut 
the  convective  time. 


4.1  Case  LarThin 

Longitudinal  profiles  of  mean  and  rms  densities  obtained  with  the  three  meshes  are  shown 
in  Fig.  6ab  for  case  LarThin.  The  fine  and  intermediate  grids  feature  results  which  are 
essentially  identical.  The  same  observations  are  made  for  radial  profiles  of  density  (Fig. 
6c),  with  additional  resolved  density  fluctuations  in  the  near  injector  region  as  the  grid  is 
refined,  as  expected  (Fig.  6d).  The  intermediate  grid  is  thus  well  suited  to  the  simulation 
of  case  LarThin  and  is  used  in  this  study. 

4.2  Case  SarThick 

The  same  study  as  for  case  LarThin  is  carried  out  for  the  case  SarThick.  Longitudinal 
profiles  of  mean  and  rms  densities  are  shown  in  Fig.  7ab.  Again,  the  intermediate  and 
fine  grid  profiles  are  close,  but  a  small  departure  is  noticeable  at  x  =  2d,,  after  the 
sudden  density  decrease.  This  region  is  characterized  by  a  strong  back  flow  and  high 
turbulence  level.  This  departure  at  x  =  2d*  is  also  present  in  the  radial  profiles  shown  in 
Fig.  7cd,  while  the  other  radial  profiles  are  similar  for  intermediate  and  fine  grids.  Since 
the  physical  phenomena  are  captured  by  both  grids  (sudden  decrease  of  density  due  to 
the  presence  of  a  back  flow)  and  differences  are  limited  and  considering  the  large  amount 
of  CPU  time  required  by  the  fine  grid,  it  is  reasonable  to  use  the  intermediate  mesh  used 
for  all  further  simulations. 
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Figure  6:  Case  LarThin  -  grid  convergence  study.  Longitudinal  profiles  of  (a)  mean, 

(b)  rms  densities  and  radial  profiles  of  (c)  mean  and  (d)  rms  densities. - coarse  — 

intermediate  □  fine 


Figure  7:  Case  SarThick  -  grid  convergence  study.  Longitudinal  profiles  of  (a)  mean, 

(b)  rms  densities  and  radial  profiles  of  (c)  mean  and  (d)  rms  densities. - coarse  — 

intermediate  □  fine 
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5  Results  without  modulation 


This  section  is  concerned  with  the  modulation  free  cases.  Results  from  simulations  are 
first  compared  with  the  available  experimental  data  in  Sec.  5.1.  and  discussed  in  Sec. 
5.2  and  5.3. 

5.1  Comparisons  with  experimental  data 

Simulations  are  first  compared  with  shadowgraph  images.  In  the  dark  core  light  is  devi¬ 
ated  by  changes  in  the  refraction  index.  This  in  turn  delineates  regions  where  the  density 
takes  large  values.  It  is  generally  admitted  that  the  dark  region  represents  the  position 
of  the  dense  core  and  highlights  its  external  boundary  giving  an  insight  into  the  main 
structures  present  in  the  flow.  Because  the  shadowgraph  images  integrate  deviations  over 
the  line  of  sight,  it  is  natural  to  integrate  density  or  density  gradients  obtained  numeri¬ 
cally  over  the  line  of  sight  and  compare  these  integrated  maps  with  experimental  data. 
Whatever  the  method,  the  comparison  will  remain  essentially  qualitative. 

Results  for  case  LarThin  are  shown  in  Fig.  8.  The  simulation  properly  retrieves  the  jet 
spreading  rate  and  is  in  good  qualitative  agreement  in  terms  of  jet  penetration. 

Case  SarThick  is  shown  in  Fig.  9.  Simulations  retrieve  a  short  “dark  core” ,  as  observed 
experimentally.  Moreover,  the  density-integrated  visualization  properly  reproduces  the 
jet  pinching  observed  experimentally.  This  comparison  shows  that  the  main  feature  of 
the  flow  are  well  retrieved  by  the  simulation.  Since  density-integrated  results  are  in  good 
agreement  with  experimental  images  for  both  cases  Larthin  and  SarThick,  they  will  be 
used  for  all  other  experimental  comparisons. 

Finally,  radial  profiles  of  average  temperature  at  one  inner  injector  diameter  are  com¬ 
pared  with  experimental  measurements  in  Fig.  10.  Simulation  results  properly  follow 
experimental  curves  for  case  LarThin:  the  plateaux  are  well  retrieved  and  slopes  are  in 
good  agreement.  Case  SarThick  shows  reasonable  agreement  with  experiment,  but  there 
are  some  small  departures  in  the  outer  region.  They  cannot  be  attributed  to  a  limited 
grid  resolution,  as  shown  by  the  grid  convergence  study.  At  this  point  in  time,  the  reason 
for  this  small  departure  is  not  found,  but  may  be  due  to  the  adiabatic  wall  condition 
used  at  the  lip.  However,  the  study  of  heat  transfer  at  the  lip  is  out  of  the  scope  of  the 
present  work. 
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(a) 


(c) 

Figure  8:  Case  LarThin  Qualitative  comparison  between  backlighting  images  from  (b) 
experiment  and  (a)  (c)  simulation.  Density  is  integrated  along  the  line  of  sight  in  (a), 
while  density  gradient  is  integrated  in  (c). 
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(c) 

Figure  9:  Case  SarThick  Qualitative  comparison  between  backlighting  images  from  (b) 
experiment  and  (a)  (c)  simulation.  Density  is  integrated  along  the  line  of  sight  in  (a), 
while  density  gradient  is  integrated  in  (c). 


Figure  10:  Radial  profile  of  temperature  at  one  injector  diameter  (a)  Case  LarThin  (b) 
Case  SarThick. 
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5.2  Instantaneous  flow  visualizations 

Longitudinal  slices  of  axial  velocity  are  shown  in  Fig.  11  for  cases  LarThin  and  SarThick. 
The  outer  jet  is  injected  at  a  higher  velocity  than  the  inner  jet,  producing  a  shear  between 
the  two  streams  that  promotes  the  inner  jet  destabilization.  The  inner  high  density  jet  is 
then  rapidly  mixed  with  the  surrounding  gas.  In  both  cases,  the  flow  is  fully  turbulent. 
However,  the  flow  structure  differs  between  the  two  cases.  While  for  case  LarThin,  the 
outer  jet  qualitatively  follows  the  inner  jet  centerline  axis  direction,  it  is  quickly  bent 
toward  the  inner  jet  injection  plane  for  case  SarThick.  In  this  latter  case,  a  recirculation 
region  is  present  behind  the  high  density  jet,  which  explains  the  sudden  density  drop 
observed  in  Fig.  7.  A  large  scale  recirculation  is  established  behind  the  inner  thick  lips, 
where  mass  entrained  from  the  inner  jet  mixes  with  hotter  gas. 


Figure  11:  Longitudinal  slice  of  axial  velocity  with  iso-lines  of  density  (100,  300  and  500 
kg  m-3)  (a)  Case  LarThin  (b)  Case  SarThick. 


Instantaneous  iso-surfaces  of  density  are  shown  in  Fig.  12  and  13  for  cases  LarThin  and 
SarThick.  Turbulence  strongly  wrinkles  the  jet  surface.  Fig.  lib  and  12b  qualitatively 
show  less  wrinkling  than  the  one  observed  on  Fig.  11a  and  12a,  suggesting  that  the  lower 
density  regions  are  more  affected  than  the  higher  density  part  of  the  flow.  Fig.  12b  and 
13b  highlight  the  major  differences  in  terms  of  jet  geometry  between  the  two  cases.  The 
flow  is  confined  by  the  back  flow  region  surrounding  the  inner  jet  for  case  SarThick  and 
cold  nitrogen  recirculates  behind  the  thick  lip  (Fig.  13b). 

5.3  Average  fields 

Longitudinal  slices  of  mean  and  rms  densities,  mean  and  rms  axial  velocities  are  shown 
in  Figs.  14  and  16  for  case  LarThin  and  in  Figs.  15  and  17  for  case  SarThick.  For  both 
cases,  a  recirculation  is  generated  behind  the  inner  and  outer  lips.  This  region  completely 
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(a) 


(b) 
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Figure  12:  Case  LarThin.  Density  iso-surfaces  colored  by  axial  velocity  (a)  86  kg  m 
(b)  300  kg  m-3. 


Axial  Velocity  [m/s] 

2.0 

0.8 

-0.5  I 

-1.8 

-3.0  U 

(b)  -■ 

Figure  13:  Case  SarThick.  Density  iso-surfaces  colored  by  axial  velocity  (a)  86  kg  m 
(b)  400  kg  m-3. 
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surrounds  the  inner  jet  in  the  SarThick  case.  As  a  consequence,  cold  nitrogen  is  entrained 
from  the  inner  jet  and  deviated  towards  the  recirculation  region  behind  the  lips  in  the 
SarThick  case,  while  the  inner  stream  propagates  downstream  in  LarThin  case.  Plots  of 
axial  velocity  are  shown  in  Fig.  18,  where  the  recalcirculation  regions  are  clearly  apparent 
for  case  SarThick  in  contrast  with  the  case  LarThin. 


Figure  14:  Case  LarThin.  Longitudinal  slices  of  (a)  mean  axial  velocity  and  (b)  rms  axial 
velocity. 
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Figure  15:  Case  SarThick.  Longitudinal  slices  of  (a)  mean  axial  velocity  and  (b)  rms 
axial  velocity. 
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Figure  16:  Case  LarThin.  Longitudinal  slices  of  (a)  mean  density  and  (b)  rms  density. 
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(b)  n 


Figure  17:  Case  SarThick.  Longitudinal  slices  of  (a)  mean  density  and  (b)  rms  density. 


Figure  18:  Longitudinal  profiles  of  (a)  mean  axial  velocity  and  (b)  rms  axial  velocity, 
Radial  profiles  at  1  inner  injector  diameter,  (c)  mean  axial  velocity  and  (d)  rms  axial 
velocity.  —  Case  LarThin - case  SarThick. 
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6  Modulated  cases 


The  modulation  free  cases  described  previously  were  essentially  meant  to  further  validate 
the  numerical  simulation  methodology.  It  is  now  possible  to  examine  modulated  cases. 
This  section  begins  with  a  description  of  the  modulation  procedure.  Results  of  simulations 
are  then  compared  with  experimental  data. 

6.1  Modulation  procedure  and  simulated  cases 

An  acoustic  modulation  is  obtained  by  imposing  harmonic  modulations  of  the  normal  ve¬ 
locity  at  the  outer  boundaries  using  a  method  devised  by  Rey  et  al.  [13] .  The  modulation 
has  the  form  : 


v\  =  Asin(27r/f)  V2  =  Asin(27r/t  +  7r)  (9) 

where  v\  and  rq  are  the  normal  velocities  at  the  lateral  boundaries,  A  the  modulation 
amplitude  and  /  is  the  modulation  frequency  set  equal  to  3.41  kHz.  These  boundaries 
are  essentially  acoustically  non-reflecting  (with  a  cut-off  frequency  of  100  Hz)  in  order  to 
avoid  the  emergence  of  transverse  eigenmodes.  All  other  boundaries  are  treated  as  in  the 
previous  cases  (in  the  absence  of  modulation) .  The  normal  velocities  at  the  boundaries  are 
in  phase  opposition,  the  jet  is  then  placed  in  a  velocity  anti-node  (maximum  transverse 
velocity  oscillation)  and  pressure  node  (no  acoustic  pressure  oscillations  along  the  jet 
centerline) . 

Modulated  cases  are  gathered  in  Tab.  6.  Two  modulation  amplitudes  are  considered:  1 
and  2  m  s-1.  These  amplitudes  are  higher  than  those  used  in  the  experiment  and  are 
chosen  to  highlight  the  driving  mechanisms. 


Case 

Geometry  and  injection  conditions 

f  [Hz] 

u’ac  [m  s-1] 

LarThin2 

LarThin 

3  410 

2 

LarThin1 

LarThin 

3  410 

1 

SarThick2 

SarThick 

3  410 

2 

SarThick1 

SarThick 

3  410 

1 

Table  6:  Injection  characteristics  of  the  modulated  cases.  u'ac  =  2  A  is  the  acoustic  velocity 
modulation  amplitude  along  the  jet  centerline. 


6.2  Comparisons  with  experimental  data 

Instantaneous  results  are  compared  with  experimental  data  for  case  LarThin  in  Fig.  19 
by  integrating  the  density  over  the  line  of  sight.  One  notices  that  the  sinusoidal  motion 
observed  in  the  experiment  is  well  retrieved  by  the  simulations. 

Comparisons  for  case  SarThick  are  shown  in  Fig.  20.  Again,  the  jet  behavior  is  qualita¬ 
tively  recovered.  Both  experiment  and  simulation  feature  an  abrupt  termination  and  t 
the  jet  pinching  observed  previously  without  acoustic  modulation  is  no  longer  present. 
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Figure  19:  Case  LarThin 1  Qualitative  comparison  between  backlighting  images  from  (b) 
experiment  and  simulation  (a).  Numerical  results  show  density  integrated  over  the  line 
of  sight. 


6.3  Averaged  fields 

6.3.1  Case  LarThin 

Longitudinal  slices  showing  the  mean  and  rms  densities  are  given  in  Fig.  21.  The  jet 
spreading  rate  is  strongly  increased  compared  to  the  unmodulated  case  (Fig.  16)  and 
density  fluctuations  are  now  very  high  in  the  near  vicinity  of  the  injector  exhaust.  The 
axial  velocity  (Fig.  22)  is  quickly  disturbed  by  the  acoustic  modulation,  and  strong 
velocity  fluctuations  are  created  near  the  injector  exit  due  to  the  formation  of  large  scale 
vortices  which  are  alternatively  released.  Both  in  terms  of  density  and  axial  velocity,  the 
flow  is  more  rapidly  homogenized  than  in  the  unmodulated  case,  indicating  that  mixing 
is  improved. 

The  acoustic  modulation  induces  an  important  reduction  of  the  jet  penetration,  which 
can  be  quantified  from  Fig.  23.  The  jet  intact  core  length  (distance  from  the  injector 
of  the  section  at  which  the  density  reaches  99  %  of  the  injection  density)  is  reduced  to 
0.30  of  the  unmodulated  case  for  case  LarThin2.  Intact  core  length  reduction  falls  to  0.4 
times  the  unmodulated  length  for  case  LarThin1.  The  principal  characteristic  values  are 
gathered  in  Tab.  7. 
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Figure  20:  Case  SarThick 1  Qualitative  comparison  between  backlighting  images  from  (b) 
experiment  and  simulation  (a).  Numerical  results  show  density  integrated  over  the  line 
of  sight. 
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Figure  21:  Case  LarThin2.  Longitudinal  slice  of  (a)  mean  density  and  (b)  rms  density. 
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Figure  22:  Case  LarThin 2 .  Longitudinal  slice  of  (a)  mean  axial  velocity  and  (b)  rms  axial 
velocity. 


Figure  23:  Longitudinal  profiles  of  (a)  mean  density,  (b)  rms  density,  (c)  mean  axial 

velocity  and  (d)  rms  axial  velocity.  —  Case  LarThin2  •  •  •  Case  LarThin1 - Case 

LarThin 
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Case 

LarThin 

LarThin1 

LarThin2 

T  ic 

^ p 

T  ic  /  T  ic, 0 
LP  /  LP 

to 

OO 

a- 

0.98  d i 
0.40 

0.73  d, 
0.30 

Table  7:  Intact  core  length  L);  for  cases  LarThin,  LarThin2  and  LarThin1.  L“:,°  is  the 
intact  core  length  for  the  unmodulated  case. 

6.3.2  Case  SarThick 

Case  SarThick  less  affected  by  the  modulation  than  the  case  LarThin.  The  density  field 
keeps  the  same  geometry  (Fig.  24)  as  that  found  in  the  modulation  free  case  (Fig.  17), 
with  a  short  high  density  core  and  the  recirculation  of  cold  nitrogen  behind  the  lips.  The 
size  of  the  recirculation  region  that  surrounds  the  inner  jet  is  nevertheless  augmented 
(Fig.  25). 


(b)  n 


Figure  24:  Case  SarThick?.  Longitudinal  slice  of  (a)  mean  density  and  (b)  rms  density. 


Due  to  the  acoustic  modulation,  the  intact  core  length  of  the  inner  jet  is  once  again 
notably  reduced  compared  with  the  unmodulated  case  (Fig.  26).  The  jet  intact  core 
length  is  reduced  to  0.40  of  the  unmodulated  case  for  case  SarThick2  and  to  0.51  for  case 
SarThick1.  Intact  core  lengths  are  gathered  in  Tab.  8. 
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Figure  25:  Case  SarThick.  Longitudinal  slice  of  (a)  mean  axial  velocity  and  (b)  rms 
axial  velocity. 


Figure  26:  Longitudinal  profiles  of  (a)  mean  density,  (b)  rms  density,  (c)  mean  axial 

velocity  and  (d)  rms  axial  velocity.  —  Case  SarThick2  •  •  •  Case  SarThick1 - Case 

SarThick 
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Case 

SarThick 

SarThick1 

SarThick2 

T  ic 

^ p 

T  ic  /  T  ic, 0 

lp  /  lp 

1.52  d, 

1 

0.77  d i 
0.51 

0.61  d, 
0.40 

Table  8:  Intact  core  length  L“:  for  cases  SarThick,  SarThick2  and  SarThick1.  Ljf’0  is  the 
intact  core  length  for  the  modulation  free  case. 

6.4  Flow  dynamics 

6.4.1  Case  LarThin 

Instantaneous  iso-surfaces  of  density  (86  kg  m-3)  over  a  half  period  are  shown  in  Fig.  27. 
The  transverse  motion  of  the  inner  jet  starts  at  the  injector  exit  and  is  amplified  further 
downstream.  The  jet  quickly  breaks  down  into  several  isolated  small  pockets  that  are 
eventually  mixed  with  the  surrounding  hotter  gas.  Vortices  are  alternatively  created  at 
the  injector  exit  (Fig.  28).  They  are  quickly  destroyed  by  turbulence  further  downstream. 
The  high  density  jet  is  initially  pinched  by  the  modulation.  The  resulting  thinned  high 
density  layer  is  wrinkled  by  turbulence  and  finally  mixes  with  the  outer  jet. 

In  order  to  highlight  periodic  phenomena  and  remove  wrinkling  due  to  turbulence,  phase 
averaged  fields  are  computed  over  20  periods  of  modulation.  This  limited  number  of 
periods  is  sufficient  to  get  rid  of  most  stochastic  fluctuations. These  fields  are  shown  in 
Fig.  29.  The  acoustic  modulation  generates  large  scale  counter  rotating  vortices  that 
stretch  the  inner  jet,  creating  the  thin  layer  observed  in  Fig.  28.  This  phenomenon  can 
be  properly  visualized  in  Fig.  33.  The  acoustic  flow  circumvents  the  high  density  stream 
and  recirculates  behind  the  high  density  inner  jet.  The  lower  density  stream  is  notably 
influenced  by  the  acoustic  forcing  and  is  peeled  away  from  the  inner  jet  forming  the  large 
scale  vortices  than  can  be  seen  in  Fig.  29. 

This  dynamical  behavior  induces  the  complex  motion  observed  in  Fig.  27.  Coupling  Large 
Eddy  Simulations  and  phase  locked  averaging  techniques  makes  it  possible  to  decompose 
this  motion  into  the  elementary  processes  that  are  observed  in  Fig.  29. 
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Figure  27:  Case  LarThin2.  Density  iso-surface  (86  kg 
period  (a)  <p  =  0  (b)  (p  =  n/2  (c)  <p  =  tx  (d)  (p  =  3n/2. 
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Figure  28:  Case  LarThin2.  Longitudinal  slices  of  instantaneous  density  during  4  phases 
(f)  of  one  period  (a)  </>  =  0  (b)  0  =  7t/2  (c)  0  =  n  (d)  </>  =  37t/2. 
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Figure  29:  Case  LarThin2.  Phase  averaged  density  iso-surfaces  (86  kg  m  3)  during  4 
phases  0  of  one  period  (a)  0  =  0  (b)  0  =  7t/2  (c)  0  =  n  (d)  0  =  3ir/2. 
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Figure  30:  Case  LarThin 2 .  Phase  averaged  transverse  slices  of  density  (blue)  with  velocity 
vectors  during  4  phases  0  of  one  period  (a)  0  =  0  (b)  0  =  7t/2  (c)  0  =  7 r  (d)  0  =  3ty/2. 
The  slices  are  positioned  at  x  =  di  and  x  =  2 di  from  the  injector  exit. 
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6.4.2  Case  SarThick 

As  outlined  previously  in  Sec.  6.3.2,  this  case  is  less  influenced  by  the  acoustic  modulation 
than  case  LarThin.  Examining  the  iso-surfaces  of  density  shown  in  Fig.  31  one  observes 
a  transverse  motion  of  the  outer  mixing  layer  close  to  the  outer  injection  channel.  This 
motion  is  hardly  visible  in  the  instantaneous  images  but  can  be  revealed  by  phase  averag¬ 
ing.  The  coherent  motion  induced  by  the  acoustic  modulation  is  highlighted  in  Fig.  32. 
Due  to  the  lip  thickness,  the  central  core  itself  is  hardly  perturbed  by  the  modulation. 
On  its  periphery,  an  alternate  release  of  vortices  is  generated  by  the  transverse  acoustic 
waves. 

This  is  confirmed  in  Fig.  33  where  two  cuts  of  density  at  dj/2  and  d,  are  presented  for 
four  different  phases  in  an  acoustic  cycle.  The  velocity  field  induced  by  the  transverse 
modulation  is  shown  to  promote  the  alternate  release  of  vortices. 
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Figure  31:  Case  SarThic}? .  Density  iso-surface  (100  kg  m  3)  during  4  phases  0  of  one 
period  (a)  0  =  0  (b)  (f)  =  tt/2  (c)  </>  =  7t  (d)  0  =  3tt/2. 
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Figure  32:  Case  SarThic}? .  Phase  averaged  longitudinal  slices  of  density  during  4  phases 
(f)  of  one  period  (a)  <p  =  0  (b)  <j)  =  7t/2  (c)  <p  =  ir  (d)  (p  =  37t/2. 
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Figure  33:  Case  Sar Thick2 .  Phase  averaged  transverse  slices  of  density  (blue)  with  veloc¬ 
ity  vectors  during  4  phases  0  of  one  period  (a)  0  =  0  (b)  0  =  7t/2  (c)  0  =  7 r  (d)  0  =  3tt/2. 
The  slices  are  positioned  at  x  =  0.5  <7,  and  x  =  di  from  the  injector  exit. 
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7  Conclusions 


Results  of  Large  Eddy  Simulations  of  two  experimental  configurations  from  AFRL,  Ed¬ 
wards,  are  reported  in  this  document.  These  two  injection  systems  operate  at  high  pres¬ 
sure  above  the  critical  value  requiring  real  gas  representations  of  the  equation  of  state 
and  thermodynamics.  The  two  experimental  conditions  feature  large  differences  in  terms 
of  flow  dynamics  and  response  to  an  acoustic  modulation.  In  the  case  of  a  large  aspect 
ratio  injector  featuring  a  thin  lip  separating  the  coaxial  streams  (LarThin),  the  flow  is 
close  to  a  classical  coaxial  jet  geometry.  In  the  second  case  (SarThick)  characterized  by 
a  small  aspect  ratio  and  a  thick  lip  between  the  two  injected  streams,  a  large  back  flow 
region  is  established  around  the  inner  jet.  This  region  strongly  confines  the  inner  jet 
and  enhances  mixing.  Simulations  are  compared  with  experimental  data  in  the  absence 
or  in  the  presence  of  acoustic  modulation.  It  is  shown  that  the  jet  geometry  is  suitably 
retrieved  both  in  absence  of  modulation  or  under  transverse  acoustic  forcing.  Simulations 
under  acoustic  modulation  capture  the  main  coupling  phenomena.  Under  acoustic  mod¬ 
ulation,  the  case  LarThin  features  a  large  scale  flag-like  motion  with  an  important  jet 
length  reduction.  The  acoustic  modulation  produces  large  scale  counter-rotating  vortices 
which  stretch  the  inner  jet  and  enhance  the  mixing  process.  In  the  case  SarThick,  the 
acoustic  modulation  has  a  limited  effect  on  the  coaxial  jet  geometry,  probably  because 
of  the  presence  of  the  dynamical  back  flow  that  dominates  mixing  in  the  injector  near 
field.  It  is  found  however  that  modulation  has  an  important  effect  on  the  jet  penetration 
length.  It  is  suggested  that  large  eddy  simulations  can  be  used  to  guide  injector  selection 
which  could  lead  to  a  lesser  sensitivity  to  external  modulation  which  could  be  beneficial 
to  a  stable  operation  of  the  propulsion  system. 
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Nomenclature 


Injector 

Si 

Se 

Si 

di 

he 
A  t 

Tc 

A x,min 


and  mesh  characteristics 

inner  lip  thickness 
outer  injector  section  surface 
inner  injector  section  surface 
inner  injector  diameter 
outer  injector  diameter 
time  step 
convective  time 

minimum  characteristic  cell  size 


Injection  conditions 


T 

inner  injector  injection  temperature 

T 

outer  injector  injection  temperature 

Ui 

inner  injector  injection  bulk  velocity 

ZLe 

outer  injector  injection  bulk  velocity 

Pi 

inner  injector  injection  density 

Pe 

outer  injector  injection  density 

J  =  {peul)/{piufj 

momentum  flux  ratio 

M  (Pe^e)/ {.Pi'U'i ) 

mass  flux  ratio 

S  Pe/  Pi 

density  ratio 

U  ZLe  /  Hi 

velocity  ratio 
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Governing  equations  and  thermodynamics 

w  conservative  variables  vector 
p  Density 

u  velocity  vector 

E  total  energy 

es  internal  energy 

G}.  kinetic  vector 

p  pressure 

x  spatial  coordinate  vector 

x  axial  distance  from  injector 

t  time 

t  viscous  stress  tensor 

q  heat  flux  vector 

r*  subgrid  scale  stress  tensor 

qf  subgrid  scale  heat  flux  vector 

A  subgrid  scale  conductivity 

Cp  heat  capacity  at  constant  pressure 

Prt  turbulent  Prandlt  number 

pt  turbulent  dynamic  viscosity 

T  temperature 

R  perfect  gas  constant 

W  molar  mass 

Tc  critical  temperature 

pc  critical  pressure 

u  acentric  factor 

Tr  reduced  temperature 


Acoustic 

/  acoustic  modulation  frequency 

u'ac  peak  to  peak  acoustic  velocity  modulation  amplitude 
<f>  phase  in  a  modulation  period 


Jet  length  scales 

Llp  intact  core  length 

Ly  ]  intact  core  length  of  the  unmodulated  case 
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